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Abstract 


‘The spectral leakage (SL) from windowing and the picket fence effect (PEF) from discretization’ have been 
among the standard contents in textbooks for many decades. The SL and PEF would cause the distortions in 
amplitude, frequency, and phase of signals, which have always been of concern, and attempts have been made to 
solve them. This paper proposes two novel decomposition theorems that can totally eliminate the SL and PEF, they 
could broaden the knowledge of signal processing. First, two generalized eigenvalue equations are constructed for 
multifrequency discrete real signals and complex signals. The two decomposition theorems are then proved. On 
these bases, exact decomposition methods for real and complex signals are proposed. For a noise-free 
multifrequency real signal with m sinusoidal components, the frequency, amplitude, and phase of each component 
can be exactly calculated by using just 4m-—1 discrete values and its second-order derivatives. For a multifrequency 
complex signal, only 2m—1 discrete values and its first-order derivatives are needed. The numerical experiments 
show that the proposed methods have very high resolution, and the sampling rate does not necessarily obey the 


Nyquist sampling theorem. With noisy signals, the proposed methods have extraordinary accuracy. 
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1. Introduction 


Since Cooley and Tukey proposed the fast Fourier transform (FFT) of Fourier analysis in 1965 [1, 
2], the technique has been indispensable in electronics, communication, signal analysis, digital 
image and audio processing, and many other fields [2, 3]. However, while performing the FFT 
algorithm, distortions in amplitude, frequency, and phase — caused by spectrum leakage (SL) due 
to signal truncation and the picket fence effect (PEF) due to frequency discretization — are 
inevitable, which have always been of concern, and attempts have been made to solve them. 

In 1970, Rife and Vincent studied the correction measurement of frequencies and levels of tones 
using discrete Fourier transform (DFT) [4]. A number of scholars have also studied this issue 
[4-18]. Interpolation techniques are the most studied and applied methods in many engineering 
fields [4-10]. Interpolation techniques refer to interpolated discrete Fourier transform (IpDFT) or 
interpolated fast Fourier transform (IpFFT). The IpDFT and IpFFT reduce SL and PEF through 
windowing and interpolation, respectively. Generally, windows with the maximum side lobe decay 
efficiently reduce SL; thus, the Hanning window is one of the commonly used [5, 7, 8]. According 
to the distribution of spectrum lines in the windows, weighting the discrete signals before DFT 
and FFT are performed, and thus the PEF is decreased to a minimum, and accurate amplitudes, 


frequencies, and phases of the signal components can be computed [4-10]. IpDFT and IpFFT can 
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be quickly and easily implemented; however, their results are not accurate when sinusoids are not 
well separated in frequency [5-8]. Weighted phase average (WPA) is another popular choice of a 
correction technique for FFT [11-15]. Relying on weighted linear regression on the phases, WPA 
averages the weighted phase estimators that are obtained by FFT within different nonoverlapping 
segments of the signals. WPA has the advantage of accurately calculated frequencies and phases, 
but amplitudes are dependent on the windows. There are many other correction techniques, such 
as the phase difference correction method [16] and sliding-window DFT method [17, 18, 27]. All 
of the correction techniques have advantages and disadvantages. 

Distinct from techniques based on DFT and FFT, matrix-based singular value decomposition 
(SVD) and singular spectral decomposition (SSD) are techniques that aim to represent signals as a 
linear superposition of elementary variable modes with unnecessary harmonic components [19— 
22]. The aforementioned techniques do not provide estimators for the spectrum, but they are 
powerful denoising filters capable of separating autocoherent features, such as anharmonic 
oscillations and quasiperiodic phenomena, from random features. They are non-parametric 
techniques. 

Other techniques commonly used to extract signal features include the wavelet transform (WT) 
[23, 24], Hilbert-Huang transform (HHT) [25], estimation of signal parameters via rotational 
invariance technique (ESPRIT) [26, 27], and multiple signal classification (MUSIC) [28, 29]. 
However, it is impossible to exactly compute frequencies, amplitudes, and phases of 
multifrequency discrete signals using the aforementioned techniques. 

This paper constructs the generalized eigenvalue equations and proves the decomposition 
theorems for multifrequency discrete real signals and complex signals. Based on the two 
decomposition theorems, decomposition methods for real and complex signals are proposed. For a 
noise-free real signal with m sinusoidal components, the frequency, amplitude, and phase of each 
component can be exactly calculated by using just 4m—1 discrete values and its second-order 
derivatives. For a complex signal, only 2m—1 discrete values and its first-order derivatives are 
needed. The decomposition results are exact in theory. The numerical experiments show that the 
proposed methods have very high resolution, and the sampling rate does not necessarily obey the 


Nyquist sampling theorem. With noisy signals, the proposed methods have extraordinary accuracy. 


2. Decomposition theorem of multifrequency discrete real signal 


Multifrequency sinusoidal signal is one of the most common signals in engineering. The 
Multifrequency real signal can be presented as follows 


m 


x(t) = > a; sin(2rf;t + Pio) (1) 


i=1 


where ai, fi, and Pio are the amplitude, frequency and phase of the ith component of signal 
x(t); m is the number of component signals. 
Discretize x(t), and construct the following generalized eigenvalue equation of the 


multifrequency discrete real sinusoidal signal 
—D,v = AXv (2) 


where A is the generalized eigenvalue; v is the generalized eigenvector; X is the square Hankel 


matrix of x(t); and D, is the square Hankel matrix of the second-order derivative X(t). And 
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where Xg, g (k = 1,-::,2n — 1) are the discrete series of x(t) and #(t), respectively. 
Theorem 1: When n > 2m, the generalized eigenvalue equation of the multifrequency real 


sinusoidal signal has 2m nonzero generalized eigenvalues, which are given by the following 


equations 
Agi-1 = Aai = 4n? ff, (C= 1m) (4) 
2.1 Proof 


The square Hankel matrix X of signal x(t) is rewritten as follows 
m 
X=Y ax; 6) 
i=1 
where a;X; is the square Hankel matrix of the ith sinusoidal component. The kth row and Ith 


column element of matrix X; is given by the following equation 


xt, = sin(2nf;i(k +l- 2)/fs + Gio), (k,l =1,--,n) (6) 
were f; is sampling rate. The square Hankel matrix Dy is given by the following equation 
m m 
D, = > Dy = ) —4n? ff aiXi (7) 
i=1 i=1 
The ranks of matrices X and D, are less than or equal to 2m, because the ranks of the 
square Hankel matrices X; and D,, are two [30]. 
Suppose A = 4r? f,7, and considering the relation Dy; = —4r? f? aiX;i, eq. (2) is modified to 


obtain the following equation 


m m 
> Dyk + An? f? 7 atao =0 (8) 
k=1,k#i k=1,k#i 


Eq. (8) has at least two nonzero solutions, which are v2; and V2;_,, because it is a 
homogeneous equation about vector v, and the rank of the coefficient matrix is less than or equal 
to 2m—2. 

In eq. (8), replacing v with v.;, and left multiplying by v3,, yields 


m m 
4r? f? = vi ( ` Dava onf ` ae | (9) 
kiki k=Lkel 


where v3, is the transpose vector of v;. Considering the following identical equation 


T 
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T 
V2, A,X (V2; 


The following equation is derived from eq. (9) and (10). 
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Eq. (11) reveals that 41?f?% is a Rayleigh quotient of matrices —D, and X, that is, 
à = 4r? f? is a generalized eigenvalue of eq. (2), and vz; is the corresponding generalized 
eigenvector [29]. Replacing v2; with Vz;-1, eq. (11) is also established. Therefore, 2 = 4r? f? 
is a double eigenvalue of eq. (2). 


Decomposition theorem 1 of the discrete real signal is thus proven. 
2.2 Component amplitudes and phases of multifrequency discrete real signal 


After the ith double eigenvalue and its two corresponding eigenvectors are computed, considering 


eq. (6), V5,a;X;V2; and v5;_,a;X;V2;-, can be expanded as follows 


(12) 
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where H;, and Hj. are square matrices with dimensions n. Their kth row and lth column 


elements are given by the following equation 


fee = sin(2mfj(k +1 —2)/f,) 


{A} = cos(2nf;(k +1 —2)/f,)’ (k,l = 1,--n) E 


For the ith generalized eigenvector, v},a;X;v2; can be replaced with v3;Xv;, and 
V5;-14;:X{V2i-1 can be replaced by Vli 1XV2i-1, because the generalized eigenvectors are 
weighted orthogonally, that is, V5; X Vk = V3;_,XV,_ = 0 (k = 2i, k + 2i — 1). Thus eq. (12) is 


modified to create 
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The amplitude a; and phase jọ can be calculated after the unknowns a;cos@jg and 
a;Sin@jg are computed from eq. . The frequencies, amplitudes, an ases of all sinusoida 
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component signals can be computed in this way. 


3. Decomposition theorem of multifrequency discrete complex signal 


The multifrequency complex signal can be expressed as follows 
m 
y(t) = > a; eJ2nfit+ io) (15) 
i=1 
where ai, fi, and Pio are the amplitude, frequency and phase of the ith component of complex 
signal y(t); m is the number of component complex signals. 
Discretize y(t), and construct the generalized eigenvalue equation of the multifrequency 


discrete complex signal as follows 
—jDyv = 2Yv (16) 


where / is the generalized eigenvalue; v is the generalized eigenvector; Y is the square Hankel 
matrix of y(t); and D,, is the square Hankel matrix of the first-order derivative y(t). And 
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where Vx, Yg (k = 1, ,2n — 1) are the discrete series of y(t) and y(t), respectively. 
Theorem 2: When n > m, the generalized eigenvalue equation of the multifrequency 
complex signal has m nonzero generalized eigenvalues, which are given by the following 


equation 
A, = 2nf,, (i= 1,-+,m) (18) 
3.1 Proof 
The square Hankel matrix Y of complex signal y(t) can be rewritten as follows 
Y= p ay; (19) 
i=1 


where a,Y; is the square Hankel matrix of the ith component complex signal. The kth row and 


lth column element of matrix Y; is given by the following equation 
yi, = eF@rfilktl-2)/fst pio), (k,l = 1,-,n) (20) 


were f; is sampling rate. The square Hankel matrix D, is given by the following equation 


M: 


m 
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The ranks of matrices Y and D, are less than or equal to m, because the ranks of the square 
Hankel matrices Y; and Dy, are both one [30]. 
Suppose A = 2rf;, and considering the relation Dy; = j2rfia;iY;, eq. (16) is modified to 


obtain the following equation 


m 


k=1,kK#i k=1,k#i 


m 


oats =0 (22) 


Eq. (22) has at least one nonzero solution, which is v;, because it is a homogeneous equation 
about vector v, and the rank of the coefficient matrix is less than or equal to m — 1. 


In eq. (22), replacing v with v;, and left multiplying by v7, yields 


2nf, = —v) ( » Dx} farf >: cute) (23) 


k=1,k4i k=1,k#4i 


where v? is the transpose vector of v;. Considering the following identical equation 


aT 
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The following equation is derived from eq. (23) and (24). 
—yT iD,v; 
2nf, = Pi] Yi (25) 


T 
vi Yvi 


Eq. (25) reveals that 27f; is a Rayleigh quotient of matrices —jDy and Y, that is, A = 27fj 
is a generalized eigenvalue of eq. (16), and v; is the corresponding generalized eigenvector [29]. 
Therefore, A = 27f; is a eigenvalue of eq. (16). 


Decomposition theorem 2 of the discrete complex signal is thus proven. 
3.2 Component amplitudes and phases of multifrequency discrete complex signal 


After the ith eigenvalue and its corresponding eigenvector are computed, considering eq. (20), 


v} a;¥;v; can be expressed as follows 
v'a;Y ivi = ajel?iny! Giv (26) 


where G; is a square matrix with dimensions n. Its kth row and lth column element is given by 


the following equation 
igh, = ef nf i(k+l-2)/fs) (k, l= 1, AP n) (27) 


For the ith generalized eigenvector, vj a;¥;v; can be replaced with v/Yv;, because the 
generalized eigenvectors are weighted orthogonally, that is, v7) Yu, = 0 (k + i). Thus eq. (26) is 


modified to create 


ajes?i0 = i r (i = 1,---m) (28) 


Eq. (28) is a complex equation, so the amplitude a; and phase gjg (i = 1,::-m) can be 


calculated. 
4. Numerical experiments 


4.1 Decomposition of noise-free multifrequency real signal and complex signal 


Based on eq. (1) and eq. (15) as the noise-free multifrequency real signal and complex signal, 
respectively, the known frequencies, amplitudes and phases of all component signals for m=10 are 
listed in Table 1. Since the highest frequency of the component signal is 290 Hz, we can set the 
sampling rate f;=299 Hz, which is slightly higher than the highest frequency. 


Table 1. Parameters of the known signal 


fjHz) 35 60 97 97.5 120 135 160 230 270 290 
aj 1.5 3.5 2 0.1 1.2 0.8 2.5 0.8 1 0.3 
Yjo(°) 30 50 170 230 90 145 0 330 280 0 


The minimum frequency difference between two neighboring components is only 0.5 Hz, 
the relative difference is about 0.5%, and the corresponding amplitude difference is up to 


20 times as large. 


For the real signal where n=4m—1=39, solve the generalized eigenvalue eq. (2), while ten 
double eigenvalues and twenty corresponding eigenvectors were computed. Ten component 
frequencies were calculated by eq. (4), and the component amplitudes and phases were computed 
using eq. (14). 

Fig. 1 shows the known multifrequency real signal and its ten decomposed component signals. 


The 39 discrete points on the original multifrequency signal are identified by symbol ‘o’. 


Although the minimum frequency difference between two neighboring components is only 0.5 Hz, 
the relative difference is approximately 0.5%, and the corresponding amplitude difference is up to 
20 times as large. The two neighboring component signals can be decomposed accurately. The 
sampling rate is less than the Nyquist sampling rate, which is 580 Hz — twice the highest 


component frequency. 


multifrequency real signal 


0 0.01 0.02 0.03 0.04 0.05 0.06 
time/s 


Fig. 1. Multifrequency real signal and its ten decomposed component signals. The 39 discrete points on the 
original multifrequency signal are identified by symbol ‘o’. The maximum absolute errors in frequency, amplitude 


and phase are 1.99e-09, 7.79e-10 and -4.39e-08, respectively. 


For the complex signal where n=2m—1=19, solve the generalized eigenvalue eq. (16), while ten 
eigenvalues and ten corresponding eigenvectors were computed. Ten component frequencies were 
calculated by eq. (18), and the component amplitudes and phases were computed using eq. (28). 

The decomposed results are similar to those of the real signal. With the multifrequency complex 
signal, however, the decomposed results have both real component signals and imaginary 
component signals. The discrete points on the original signal are much less than the real signal. 

The absolute computational errors in the component frequencies, amplitudes and phases are 
listed in Table 2. With the real signal, the maximum absolute errors in frequency, amplitude and 
phase are 1.99e-09, 7.79e-10 and -4.39e-08, respectively. With the complex signal, these errors are 
6.41e-07, -2.5le-07 and -6.72e-06, respectively. It should be noted that all maximum errors 
resulted from the two neighboring components with minimum frequency difference and maximum 
amplitude difference, and the other errors were much less. Therefore, these errors were considered 
to be caused by the computer rather than the method. 

Therefore, the proposed two decomposition theorems are exactly in theory. The corresponding 
decomposition methods are highly accurate and possess high resolution, and it is unnecessary for 


the sampling rate to obey the sampling theorem. In fact, it is impossible to extract accurate 


frequencies, amplitudes and phases from multifrequency signals by using other methods with so 


little discrete data. 


Table 2. The absolute computation errors 


in the component frequencies, amplitudes and phases 


Errors of real signal Errors of complex signal 

Af; (Hz) Aaj A@jo(°) Af; (Hz) Aaj; AQ jo(°) 
-5.76e-12 2.07e-12 1.22e-10 -1.99e-13 -2.00e-14 1.80e-12 
3.41e-13 2.19e-13 -7.68e-12 0 0 -3.98e-13 
9.54e-11 7.79e-10 -2.23e-09 -5.18e-09 -4.23e-08 6.92e-08 
1.99e-09 -7.77e-10 -4.39e-08 6.41e-07 -2.51e-07 -6.72e-06 
1.14e-13 4.26e-14 -4.49e-12 3.50e-11 -6.50e-13 -4.22e-10 
1.25e-12 4.68e-13 -4.56e-11 2.20e-11 -3.09e-12 -2.09e-10 
-2.56e-13 9.70e-13 -3.36e-12 9.95e-13 -1.90e-13 -6.20e-12 
-8.81e-13 -8.34e-14 2.02e-11 0 -2.00e-15 0 

3.41e-13 -9.43e-14 -8.92e-12 0 -1.30e-14 1.99e-12 
3.41e-13 -4.11e-15 -7.67e-12 0 -6.99e-15 -4.23e-12 


The color bold numbers are the maximum absolute errors of the component parameters. All 
maximum errors resulted from the two neighboring components with minimum frequency 


difference and maximum amplitude difference, and the other errors are greatly less. 


Generally, the first-order and second-order derivatives of signals are unknown. In this case, they 
can be easily obtained by the differential circuits or numerically calculated by using the discrete 


values of the original signal. 
4.2 Simulation to FMCW radar measurement 


Frequency modulation continuous wave (FMCW) radar is widely used to measure the distance of 
targets. The transmitting antenna of the radar system transmits the FMCW radio signal, the 
receiving antenna receives the reflected signal from the target, and then mixing the transmission 


signal and reflection signal yields [31] 
z(t) = AeI@rfot+n) + w(t) (29) 


where A is the received signal power and w(t) is the system noise. The frequency fp and phase 


@p are given by the following equations 


_ 2BR 30 
fy =F (30) 

_ Anf-R n AmB R? si 
Pp = c c?T, ( ) 


where f, is the chirp start frequency; B is the bandwidth of the chirp; Tẹ is the duration of the 
chirp; R is the distance to the target; c is the speed of light. 

According to eqs. (23) and (24), the distance to the target can be calculated by measuring the 
frequency fp or the phase øp. In general, the distance measurement accuracy of the phase 
method is higher. However, because of the periodic ambiguity of the phase, when the measured 
distance is long, it is difficult to use the phase method [31]. Here, we simulate to measure the 


distance of the target by the frequency method and investigate the influence of noise on the 


measurement accuracy by using the proposed method of this paper. 

Set the parameters of radar as follows: the chirp start frequency f, = 24 GHz, the bandwidth 
of the chirp B = 100 MHz, the duration of the chirp T, = 512 us, and the speed of light 
c = 3 x 108 m/s. Suppose the distance of the target R = 12.3456789 m, and the noise w(t) is 
Gaussian white noise with zero mean. With a specified signal-to-noise ratio (SNR), 1000 
simulation measurements were made, and the simulated measurement results of distance were 


shown in figures 2 to 4. 
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Fig. 2. The absolute errors and relative errors of the average distance of 1000 simulation measurements vs. SNR 


Fig. 2 shows the absolute errors and relative errors of the average distance of 1000 simulation 
measurements vs. SNR. As seen from fig. 2, even if the SNR is very low, as low as -10 dB, the 
distance errors are less than 4 mm, and the relative error is approximately 0.032%. When SNR >= 
10 dB, the absolute errors and relative errors are close to zero 

The reason for the high distance accuracy is that the simulation measured frequency, amplitude 
and phase of the signal z(t) are of high accuracy. Fig. 3 shows the comparison between the mean 
squared errors (MSEs) of the simulation measured frequency, amplitude and phase and their 
Cramer-Rao bounds (CRBs) [4]. 


10log10 ( 1/MSE ) / dB 
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Fig. 3. The MSEs and CRBs of frequency, amplitude and phase of the signal vs. SNR. The MSEs of the frequency, 
amplitude and phase are less than those of their CRBs. This indicates that the method has achieved extraordinary 
accuracy. 


In general theory, it is impossible that the MSEs of the frequency, amplitude and phase are less 
than their CRBs. The simulation results, however, show that the MSEs are indeed less than their 
CRBs. It is very difficult to explain these results at this time. On the hand, this indicates that the 
method has achieved extraordinary accuracy. 

To investigate the robustness of the method, the MES of the simulation measured distance is 
calculated by the following equation 

1N 
a= nae (Ri — R)? (32) 
i=1 
where N = 1000, is the simulation times; R; is the simulation measured distance at a certain SNR; 
R is the real distance of target. 
Fig. 4 gives the curve of o2 vs. SNR. It shows that o2 decreases rapidly with the increase of 


SNR. So the method has well robustness. 


MES of distance / mm? 
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Fig. 4. The MSE of distance vs. SNR 


5. Conclusions 


The key point of this paper is the construction of generalized eigenvalue equations for real 
signals and complex signals. For the multifrequency real signal, the generalized eigenvalue 
equation is constructed by using the square Hankel matrixes of discrete values and their 
second-order derivatives of signal. For the multifrequency complex signal, the generalized 
eigenvalue equation is constructed by using the square Hankel matrixes of discrete values and 
their first-order derivatives of signal. The aforementioned Hankel matrix can be replaced by a 
Toeplitz matrix or other matrix as long as the rank of the matrix equals the rank of the 
corresponding square Hankel matrix. The first-order and second-order derivatives of signals can 
be obtained by differential circuits or numerically calculated by using the discrete values of the 
original signal. 

The proposed two decomposition theorems are exactly in theory. For a noise-free real signal 
with m components, the proposed method can exactly compute the frequency, amplitude and 
phase of each component using only 2m-1 discrete values and corresponding second-order 
derivatives. For a noise-free complex signal, the number of discrete values decreases to 2m-1, and 
the second-order derivatives are replaced by the first-order derivatives. The proposed methods 
have very high resolution. The sampling rate is unnecessary to obey the sampling theorem. 

With noisy signals, the proposed methods have extraordinary accuracy. It is possible that the 
MSEs of the frequency, amplitude and phase are less than their CRBs. 


For the multifrequency discrete real signal and complex signal, it may be the first time to realize 


the exact decomposition of frequencies, amplitudes and phases of component signals in theory. 
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